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Abstract 


The aim of this article is to present a new method based on Lucas poly- 
nomials and residual error function for a numerical solution of fractional 
Bagley—Torvik equations. Here, the approximate solution is expanded as 
a linear combination of Lucas polynomials, and by using the collocation 
method, the original problem is reduced to a system of linear equations. 
So, the approximate solution to the problem could be found by solving this 
system. Then, by using the residual error function and approximating the 
error function by utilizing the same approach, we achieve more accurate 
results. In addition, the convergence analysis of the method is investi- 
gated. Numerical examples demonstrate the validity and applicability of 
the method. 
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1 Introduction 


Fractional differential equations have important rules in many fields of sci- 
ence and engineering. For example, in viscoelasticity [4, 3], economic growth 
model and finance [5, 16], biology [24], control theory [8, 14, 20, 21], dynam- 
ics of particle [29], electrical circuits [8], and vibration [25], some issues can 
be modeled as fractional differential equations. 
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The fractional Bagley—Torvik equation was originally introduced in 1983 
to describe the motion of an immersed plate in a Newtonian fluid [30] as 


d? A d 


3 

= U(e) + cU (2) = 0, 
dx 
where m, and A are the mass and area of the plate, respectively, r is the 
fluid density, c is the spring of stiffness, and 7 is viscosity. The solution of 
the Bagley—Torvik equation has been studied by researchers for the past two 
decades. In [26] authors applied the Adomian decomposition method for the 
solution of Bagley—Torvik equation. El-Gamel and Abd-El-Hadi [9] presented 
the Legendre-collocation method to approximate the solution of fractional 
Bagley—Torvik equations. Zolfaghari et al. [34] studied an application of the 
enhanced homotopy perturbation method to find the approximate solution of 
Bagley—Torvik equation. In [32], an integral transform method is considered 
for solving Bagley—Torvik equation. Srivastava, Shah, and Abass [28] pro- 
posed a numerical method for studying Bagley—Torvik equations based on the 
Gegenbauer wavelet together with block pulse function. In [10], authors pre- 
sented Chelyshkov—Tau as an effective tool for solving Bagley—Torvik equa- 
tion. Cenesiz, Keskin, and Kurnaz [6] solved Bagley-Torvik equations by 
using the generalized Taylor collocation method. Authors of [15] utilized hy- 
brid functions approximation, which consists of the block pulse function and 
Bernoulli polynomials, for the numerical solution of Bagley—Torvik equations. 
Zahra and Van Daele [33] used a discrete spline function and nonstandard 
Grunwald-Letnikov and weighted and shifted Grunwald—Letnikov difference 
operators to propose the solution to Bagley—Torvik equations. El-Gamel and 
Abd-El-Hadi [9], by using Legendre basis functions, reduced Bagley—Torvik 
equation to a system of linear equations and by solving this system presented 
a numerical solution to the Bagley—Torvik equation. Authors of [27] intro- 
duced the numerical solution of Bagley—Torvik based on reproducing kernel 
Hilbert space. In [31], generalized Bessel functions of the first kind are applied 
for the numerical solution of the fractional Bagley—Torvik equation. 


The outlines of the article are as follows: In section 2, we briefly introduce 
the Caputo fractional derivative, Fibonacci, and Lucas polynomials and de- 
scribe their properties. In section 3, we construct a numerical method for a 
solution of fractional Bagley—Torvik equations using Lucas polynomials and 
residual error function. In section 4, the convergence analysis of the pro- 
posed method is studied. The numerical results for some problems are given 
in section 5, and at the end, we have a brief conclusion. 


2 Basic definitions and requirements 


Definition 1. If a > 0, then the Caputo fractional derivative operator of 
order a is defined as 
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1 
T(m—a 


D° f(z) = ) fe pe 0) aid. 


wherem—1l1<a<m. 


The Caputo derivative has linear property and 


D°(c) =0, cis a constant. 


k =5, 4 
peat) =( recarye * RZ lal, REN, 
0 ifk <a], kE No, 
where No = {0,1,2,...} and N is the set of natural numbers. 


The Fibonacci polynomials F(x) and and Lucas polynomials L,,(x) are 
defined by recursive relations as 


and 


Lo(z) =2, I(x) =2, 
Ly (a) = eLy-1(a) + In-2(x), n> 2, 
respectively. Here we remark that Fibonacci and Lucas polynomials are the 


special case of Chebyshev polynomials (see [22]). Lucas polynomials have 
explicit form as 


where |x| is the largest integer less or equal to x. According to [18], the 
first derivative of Lucas polynomials can be evaluated using the Fibonacci 
polynomials as 

L! (x) = nF, (2). (1) 


Continuing this approach by repeating derivation on both sides of (1) gives 
1 (a) =nFE-V(n), > 2. 


If u(a) is a continuous function, then we can approximate u(x) by a linear 
combination of Lucas polynomials as 


u(x) & $7 ajL;(x) =L(a)A, 
j=0 
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where L = [Lo(x), Li(z),...,2m(x)| and A = [ao,a1,.-.,@m]7. Moreover, 
for the kth derivation of u(x), we have the approximation 
k 
u(x) = Sab (a). 
j=0 

Therefore, as mentioned in [11], the approximation for the kth (k > 2) deriva- 
tion of u(a) can be formulated as 

ul) (x) = nF(2)D*-1A, 
where 


F= [Fo(x), Fi(a), see »Fim(2x)], 


Dickie tt) = : d ) 


and d is an m x m matrix, which is defined as 


a. afi sn 23% iff >, 
al 0 if 7 <i. 
Further details about Lucas polynomials and application of Lucas polynomi- 


als for solving problems arising in engineering, such as ordinary and partial 
differential equations, can be found in [1, 2, 7, 12, 11, 13, 18, 19). 


3 Construction of method 


Consider the fractional Bagley—Torvik equation 
AD? f(z) + BD? f(«)+C f(z) =g(x), x € [0,1], (2) 
with initial conditions 
fO=fo, FO =f, 
or boundary conditions 
fO=fo, fM=A. 


By using the Caputo fractional derivation, (2) can be rewritten as 


AD? f(x) + va [eo rOa+C He) =9@). 
0 


T(5 
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Let the approximate estimation for the solution of (2) have the following 


form: 
M 
=~ S” ay L; (2). (4) 
j=0 


Now, by collocating at the nodes {x; : i = 1,...,M—1}, where 0 < 21 < 
-<2y_-1 <1, and utilizing (3), we get 
M 
ADastilod+ sey Loos f. O-O AOA +OL oy bsla) = ola) 


j=0 


Also, initial and boundary conditions lead to 


Sos = fo, Sato = fo, (6) 


and 
Ss = fo, Sos =f, (7) 


respectively. Hence, the combination of (5) together with (6) or (7) gives a 
system of linear equations as 


UA =b, 
where, for initial conditions, 
b = [9(21),.-.,9(tae-1); fo, fol’ 
XL; mae 
AL4_,(2i) + Ta) he (a; = t) 2 Li 4 (t)de +CL;~1(x:) 
U..= ifl<i<M-1, 
oy L;-1(0) ifi= M 
Di (0) ifi=M +1, 
and for boundary conditions, 
b = [9(a1),---,9(@m_-1), fo, fil’; 
ALG ave: By Io’ (a — t) ray 1(t)d: +CL;~1(ax;) 
— ifl<i<M-l1, 
” L;-1(0) if i= M, 
L;-1(1) ifi=M+1. 


For example, if M = 2, by using Chebyshev—Gauss—Lobatto nodes, U has 
the following structure: 
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Cc 5 _B Cc 
20 FS 2A+2V27 49% 
Ur=|2 0 2 
0 1 0 
and . . ‘ 
20 F 2A+2V2- +95 
Up=|2 0 2 
2 1 3 


Therefore, by solving the obtained system of linear equations, the approxi- 
mate solution of the fractional Bagley—Torvik equation is determined. Here, 
we present a more accurate method using the residual error function [7, 17] 
for the solution of the fractional Bagley—Torvik equation. If we display the 
error of approximation (4), as 


M 
e(x) = f(z) — > aj L; (2), 
j=0 


then the error function satisfies the fractional differential equation 


A D*e(x) + B D?e(x) + C e(x) = R(x), (8) 
where 
M B M x 

R(x) =9(2)- A) a, L7(a)-— <4 Qj x—t) 2Lii(t)d (9 
He) AY OTA Le f yAL()d (9) 

M 
= Cy ot y(e). (10) 

j=0 


The above fractional differential equation is accompanied with initial condi- 


tions 
e(0) = e’(0) = 0, (11) 


or the boundary conditions 
e(0) = e(1) = 0. (12) 


Now, we propose the approximate solution to (8)—(12) using Lucas polyno- 
mials as 


N 
e(x) = S° By L;(a), N>M. 
j=0 


By using the idea described above, we can get the an approximation for error 
function e(z). So, we obtain a better approximation 
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M N 
F(a) * So ag L;(x) + S5 BL ;(2) 
j=0 j=0 


for the numerical solution of the fractional Bagley—Torvik equation. 


4 Convergence analysis 


In this section, we argue about the convergence of the proposed method. For 
this aim, first, some requirements are given. 


Lemma 1. [1] Assume that f() is an infinitely differentiable at « = 0. Then 
f(x) can be represented by using Lucas polynomials as 


SS (yap 


where 


Lemma 2. [1] For every i > 0, the Lucas polynomials can be bounded as 
ILi(x)| < 20°, 
where go is the golden ratio. 


Theorem 1. [1] Let f(x) be defined on [0,1], and there is a positive con- 
stant A such that |f@(0)| < A’, i > 0. Moreover, suppose that f(x) has a 
representation 


f(z) = dD ckil2), (13) 


and that e(@) = 337° y,4, cLi(x) is an error of approximation function f(a) 
by Lucas polynomials of degree M. Then 

A’ cosh(2A) 

as 


and the series (13) is convergent and 


2e4° cosh(2A)(Ac)!+1 


ote le (M +1)! 


In the following theorem, we discuss the convergence of the presented 
method of the previous section. 
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Theorem 2. Let f(x) be an infinitely differentiable at x = 0, and there 
is a constant A > 0 such that |f(0)| < A‘, i > 0. If e(x) is defined as 
e(a) = f(x) — ean a,L;(a) and has a representation 


then the proposed method has the error estimation 


2e4° cosh(2A)(Aa)N+1 
(N+ 1)! 


|E(x)| < 


Proof. According to the previous section, the approximation 


M N 
i=0 1=0 


] has the error E(x) = e(x) — 7’, bi Li (a). Also, 


O(a) = fO(e ) Dade 


Since L;(x) is a polynomial of degree 7, so i 


sentation: 


(a) has the following repre- 


L® (2) — J Uo F Aj F 0 + aj,_,07~* if 7 22, 
0 if 7 <i. 


Therefore 


(ny _ J % i724, 
L; ={5 ifj <i. 


If we set P = max{|L‘(0)| : i,j = 0,...,M}, then for i = 1,...,M by 
using Theorem 1, we get 


M4; 

‘ AJ cosh(2.A)P 
|e (0)| < At + y ee 
j=0 


< (A+ cosh(2A)Pe4)’. 


Moreover, for i > M +1, we have |e“ (0)| < A‘. If we apply Theorem 1 for 
the function e(x) = f(x) — an a,L; (a), then 


2e4? cosh(2A)(Aa)N+1 
1B(2)| < eA) 
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where 
A = A+cosh(2A)Pe4. 


5 Numerical results 


In this section, some examples are presented to show the accuracy of the pro- 
posed method. These examples consist of initial and boundary conditions. 
Also, to show the accuracy and validity of the proposed method, we have a 
comparison between our approach and a number of other methods. In com- 
putations, we utilize Chebyshev—Gauss—Lobatto nodes as collocation points, 
and all of the computations have been performed in MAPLE 18 software. 


Example 1. Consider fractional Bagley—Torvik equation 


3 
D? f(x) + D? f(x) + f(a) =2° +724+14+ = 
Jn 

with the initial conditions f(0) = 1 and f’(0) = 1. This problem has the 

exact solution f(z) = 2? +a-+1. Here we take M = 6 and N = 10. We 

compare the Lucas collocation method (LCM) and Lucas collocation method 

combined with residual error function (LCM-REF) with the Chelyshkov-Tau 

method [10] and Legendre collocation method [9]. Results are given in Table 

1. Absolute errors of LCM and LCM-REF are listed in Table 2 and plotted 
in Figure 1. 


Table 1: Comparisons of the presented methods for Example 1 


x Exact solution LCM-REF LCM Chelyshkov—Tau [10] Legendre collocation [9] 


0.1 1.101000 1.101000 =: 1.101000 1.101000 1.101000 
0.25 1.265625 1.265625 = 1.265625 1.265625 1.265625 
0.5 1.625000 1.625000 = 1.625000 1.625000 1.625000 
0.75 2.171875 2.171875 = 2.171875 2.171875 2.171875 

1 3.000000 3.000000 3.000000 3.000000 3.000002 


Table 2: Absolute errors of the presented methods for Example 1 


x 0.1 0.3 0.5 0.7 0.9 


LCM 8.90000E-48  1.22000E-47 1.67000E-47 2.23000E-47 = 3.28000E-47 
LCM-REF  8.72634E-48 1.06045E-47 1.19949E-47 1.30881E-47 = 1.39519E-47 


Example 2. In this example, we study the fractional Bagley—Torvik equation 
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w 


LCM 
—— - LCM-REF 


Figure 1: The plot of absolute errors of LCM and LCM-REF for Example 1 


(48p(ax) + 7/xq(2)), 


8 3 13 uy 2 
D? f(a) + 77D? F(a) + FF) = 39050 /a 


where 


p(x) =16000a* — 32480x? + 212802? — 4746a + 189, 
q(x) =32502° — 942524 + 26488027 — 44, 


with the boundary conditions f(0) = 0, f(1) = 0. This problem has the 
exact solution 


10 25° 250 


We examine the proposed method with M = 6, N = 10. In Table 3, a 
comparison between absolute errors of Lucas collocation method combined 
with residual error function (LCM-REF), Chelyshkov-Tau method [10], Harr 
wavelets method [23], and Bessel collocation method [31] is given. In Figure 2, 
the plot of the exact solution and approximate solution, which is obtained by 
the combination of the Lucas collocation method and residual error function, 
is displayed. 


Example 3. Consider the fractional Bagley—Torvik equation 


A D? f(x) +B D? f(x) +C f(x) =9(z), 
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Table 3: Comparisons of LCM-REF for Example 2 
x  LCM-REF  Chelyshkov—Tau [10] Harr wavelets [23] Bessel collocation [31] 
0.1 2.69915E-48 5.92720E-14 6.49908E-7 1.0800E-2 
0.2 3.21281E-48 1.18400E-13 6.35657E-7 8.9595E-3 
0.3  3.72702E-48 1.77249E-13 3.71584E-7 3.7797E-3 
0.4 4.23527E-48 2.35568E-13 9.48220E-7 1.4413E-7 
0.5 4.71891E-48 2.17578E-13 1.59573E-6 1.0001E-3 
0.6 5.20793E-48 2.92504E-13 1.05494E-6 6.6150E-8 
0.7 5.69369E-48 3.82671E-13 6.34678E-7 1.2599E-3 
0.8 6.20644E-48 3.82256E-13 1.88690E-6 1.2800E-3 
0.9 6.54825E-48 2.90107E-13 3.13999E-6 2.0656E-8 


Exact solution 
© _LCM-REF solution 


Figure 2: The plot of exact and LCM-REF solutions for Example 2 


with the initial conditions f(0) = 0, f’(0) = 0. This problem has the exact 
solution 


x ae 
sa) = [Gala n)alr)dy = 5) GP Ge BO, a (-F v0) 


where G3(x) is three-term Green’s function, which is defined as 


1 (-1)* C0 B 
Gy(a) = 5 Sphere W® a (Fv), 
k=0 
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and F),,, is the Mittag-Lefller function with two parameters \ and ju, and 


(k) (f+ k)ly? 
E = k=0,1,2,.... 
Awl) Qo FROG EA moe? 


Let A=1, B= $,C =, and g(x) =8. For this case, we choose M = 30 and 
N = 40. Numerical comparisons of the proposed methods with Chelyshkov— 
Tau method [10], Legendre collocation method [9], and generalized Taylor 
collocation method [6] are listed in Table 4. In Table 5, absolute errors of 
LCM and LCM-REF are displayed. Figure 3 exhibits the comparison of 
analytical and LCM-REF solutions of this example. The plot of absolute 
errors of LCM and LCM-REF is illustrated in Figure 4. 


Table 4: Comparisons of the presented methods for Example 3 


x Exact solution LCM-REF LCM — Chelyshkov-Tau [10]  Taylor-collocation [6] Legendre collocation [9] 
0.1 0.036487 0.036486 0.036483 0.036453 0.036485 0.036471 
0.2 0.140639 0.140636 0.140632 0.140575 0.140634 0.140615 
0.3 0.307484 0.307480 0.307473 0.307403 0.307476 0.307434 
0.4 0.533284 0.533278 0.533269 0.533252 0.533271 0.533225 
0.5 0.814756 0.814749 0.814739 0.814860 0.814735 0.814661 
0.6 1.148837 1.148828 1.148816 1.149069 1.148805 1.148733 
0.7 1.532565 1.532555 1.532541 1.532870 1.532521 1.532424 
0.8 1.963029 1.963018 1.963002 1.963440 1.962974 1.962874 
0.9 2.437334 2.437322 2.437305 2.437829 2.437455 2.437134 


Table 5: Absolute errors of the presented methods for Example 3 


x 0.1 0.3 0.5 0.7 0.9 


LCM 4.23327E-6 1.16836E-5 1.82885E-5 2.41776E-5 2.93789E-5 
LCM-REF  1.78503E-6 4.92605E-6 7.71291E-6 1.01973E-5 = 1.23905E-5 
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=== Exact solution 
© _LCM-REF solution 


Figure 3: The plot of exact and LCM-REF solutions for Example 3 


0.00003 4 
0.00002 4 
LCM 
—— LCM-REF 
0.00001 4 
Oe T T T T T 
10) 0.2 0.4 0.6 0.8 1 


Figure 4: Absolute errors of LCM and LCM-REF for example 3 
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6 Conclusion 


A new numerical method using Lucas polynomials was proposed to solve the 
fractional Bagley—Torvik equation. In this approach, we expanded the exact 
solution as a finite linear combination of Lucas polynomials. Then, by us- 
ing Chebyshev—Gauss—Lobatto nodes as collocation points, the approximate 
solution was obtained. To improve the results, we applied the residual error 
function, and the error function was estimated by Lucas polynomials. So we 
can improve the results and get a more accurate approximation. Numeri- 
cal tests and comparisons with other numerical methods indicated that this 
method is reliable and has acceptable accuracy. 
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